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Abstract. The uniqueness and stability of irradiated flar- 
ing passive protostellar disks is investigated in the context 
of a simplified set of equations for the vertical height H 
as a function of radius R. It is found that the well-known 
flaring disk solution with H oc i?^/^ is not unique. Diverg- 
ing solutions and asymptotically conical (H cx R) solu- 
tions are also found. Moreover, using time-dependent lin- 
ear perturbation analysis, it is found that the flaring disk 
solution may become unstable to self-shadowing. A local 
enhancement in the vertical height alters the functional 
form of irradiation grazing angle, and causes the 'sunny 
side' of the enhancement to grow and the 'shadow side' to 
collapse in a run-away fashion. This instability operates 
in regions of the disk in which the cooling time is much 
shorter than the vertical sound crossing time, which may 
occur in the outer regions of the passive irradiated disk if 
dust and gas are sufficiently strongly thermally coupled. 
Processes that may stabilize the disk, which include active 
accretion, irradiation from above (e.g. a scattering corona) 
and low disk optical depth, are likely to operate only at 
small or at large radius. The simple analysis of this Letter 
therefore suggests that the instability may alter the fiaring 
disk structure at intermediate radii (between the actively 
accreting and fast rotating inner regions and the optically 
thin outer regions). 
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1. Introduction 

It has become widely accepted that the infrared excess 
observed in the spectral energy distribution of many T- 
Tauri stars is caused by thermal emission from a circum- 



stellar dusty disk (Adams & Shu, 1986; Adams, Lada & 
Shu, 1987). However, the spectra predicted by theoreti- 



cal models of both accretion-powered disks (Lynden-Bell 



fc Pr ingle, 1974) and flat reprocessing disks (Friedjung, 
1985| ) have an infrared spectrum that goes as XFx oc 
X-^-^^. This is steeper than the XFx oc A-°-5*°-^-° typ- 



ically obser ved fr om T-Tauri stars (Rucinski, 1985 : Ryd- 
gren & Zak, 1987 ). A natural way to explain the relatively 
warm dust at large radii (which is needed to increase the 
slope of the model spectrum), is to invoke a fiaring disk 
geometry, i n wh ich H/R cx R'' , with 7 > (Kenyon & 
Hartmann, 1987 ). Because of the angle between the disk 
surface and the star, the disk intercepts considerable stel- 
lar flux even at large radii, and therefore acquires a shal- 
lower temperature profile than the T cx R~^^^ found from 
the accreting and/or flat disk models. The flaring index 
7 can be determined sclf-consistently by equating the in- 
tercepted flux with the emitted blackbody fiux, under the 
assumption that the disk is vertically isothermal. These 
models yield 7 = 2/7 (Chiang & Goldreich, 1997, hence- 
forth CG97; D'Alessio et al. |l999| , henceforth DCHLC99). 
Apart from producing the fiatter SEDs, it was recognized 
that the optically thin surface layers of the disk will have 
much higher temperatures than the interior (Calvet et 
Malbet & Bertout |1991|; CG97). In addition to 



1991 



al 

producing an extra component in the SED, this optically 
thin layer is also able to explain certain dust features in 
emission which would otherwise be expected in absorp- 
tion. A more direct piece of evidence for the flaring nature 
of protostellar disks was pr ovided by the HST images of 
IIH30 (Burrows et al. 



19961) 



Despite the successes of the flaring disk model in ex- 
plaining observations of T-Tauri stars, there are still some 
theoretical issues that remain to be addressed. In this pa- 
per we address two questions: are the flaring disk models 
stable against perturbations in vertical height and tem- 
perature, and are they unique solutions to which a time- 
dependent evolution of the passive disk settles down? The 
issue of stability has been addressed before by DCHLC99. 
They considered perturbations in the temperature of the 
disk, and found that they damp out as they propagate in- 
ward. However, their analysis is based on the assumption 
that the disk vertical height quickly responds to changes 
in the temperature, or in other words that the vertical 
sound crossing time tdyn is much smaller than the cooling 
time ithorm- This is true for the inner regions of the disk, 
but breaks down at large radii. In this Letter we study the 
other extreme case: the case of tdyn ^ ithorm, which may 
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occur in the outer regions of the disk. We start from the 
same equations as DCHLC99 and CG97, but replace the 
static equation for the vertical disk height with a simpli- 
fied one-zone dynamic model. We derive the pertubation 
equations and dispersion relations from them, and find 
that under these circumstances the disk becomes unsta- 
ble. 

We will start the analysis in Section |^ with a descrip- 
tion of the family of solutions of a static irradiated opti- 
cally thick disk. The flared disk model mentioned above 
is just one of these solutions. The other solutions either 
diverge at, or are conical {H oc R) beyond some radius. 
In Section || we will discuss the stability under the two 
extreme conditions (ttherm < tdyn and itherm > ^dyn)- In 
Section || we discuss the consequences of this instability, 
and when it is expected to show up. 
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2. Static equations for an irradiated disk 

A disk with flaring angle a intercepts a flux of i^irr — 
a{R^/R)^ajiTf ergcm^^s^^. The flaring angle a is deflned 
as the angle between the grazing incident stellar light and 
the surface of the flaring disk, and is assumed to be small. 
If we assume the disk to be vertically isothermal, the emit- 
ted flux is Fcinit = ctrT^, where To is the temperature of 
the disk. Equating Fji-r with Tcmit yields 



1/2 



When i? the flaring angle is given by 

a = R 



dR\ R 



(1) 



(2) 



where Hg is the height of the surface of the disk above 
the midplane. For a static solution the vertical pressure 
balance equation for an isothermal disk yields a Gaussian 
density proflle of the form 



HpV2TT 



exp 



(3) 



where E is the surface density, and Hp is the pressure scale 
height 



I kT^R^ 



(4) 



where fj,g is the mean molecular weight and is the 
unit atomic mass. The surface height Hg and the pressure 
scale height Hp are related by a number of the order of 
a few, dependent on the opacity. We assume this to be 
approximately a constant x = Hs/Hp. CG97 take % ~ 4 
for their analytic calculation of the disk structure. If we 
eliminate Tc from the above set of equations we arrive at 



1 / fj-gm^GM^, 
X V kT,R, 



Hp 
R 



= R 



Hp 
dR\ R 



(5) 



Fig. 1. The family of static irradiated passive disk solu- 
tions according to Eq.(^) for C = 1 x 10^^. 

Deflning the dimensionless quantities r = R/R^ and h = 
Hp/ R, one arrives at 



dh{r) 



C 



dr r 
where C is 



h{r)\ 



X \ kT^R^ 



(6) 



(7) 



This ordinary differential equation (ODE) has a 1- 
parameter family of solutions: 



h{r) 



K 



-1/7 



(8) 



where K is an arbitrary real constant. For K < the solu- 
tions diverge at r = y^—7C/2K, for K > the solutions 
become constant beyond r = y/7C/2K, and ior K = 
one obtains a powerlaw solution. These three classes of 
solutions are shown in Fig.(|l|). 

The power law solution is given by 



h{r) 



2 

7C 



1/7 



,,2/7 



(9) 



which is the solution given by CG97. All the other so- 
lutions asymptotically tend to this power law solution, 
Eq.(^, as r J, 0. But as r goes to infinity the solutions 
deviate from this power law at some radius, and either 
diverge to /i^oo at a particular radius, or asymptotically 
approach a conical geometry /i(r)^ const, i.e. Hg cx R. 

The direction in which the solutions deviate from each 
other is interesting. Since the inner regions of the disk 
(small radii) cast their shadow over the outer regions 
(large radii), the "causality" of the system is pointing out- 
wards. A change in the structure of the disk at small radii 
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has its repercussions on the disk at large radii, but not 
vice versa. In principle this means that the (numerical) 
integration of the ODE {Eqjs]) must be done from inside- 
to-outside (i.e. in causal direction). However, in this di- 
rection the solutions tend to strongly deviate from each 
other. A slight change in integration constant at the start 
of the integration causes a very large difference at large 
radii. For stable integration of the ODE (Eq.|^) one should 
place the boundary condition at the outside and integrate 
inwards. Though mathematically correct, it is physically 
meaningless, and may hint to an intrinsic instability of the 
time-dependent equations. 



3. Time-dependent equations for irradiated disk 

The stability of the flaring (power law) solution can be 
studied analytically using time-dependent linear pertur- 
bation analysis. We study two regimes, one of which has 
been studied before by DCHLC99. 



3.1. Temperature perturbations at hydrostatic equilibrium 

In the analysis of DCIILC99 the cooling time scale of the 
disk was assumed to be much larger than the orbital time 
scale, which is true for the inner regions of the disk. We 
may then assume that the disk will always be in hydro- 
static equilibrium, and we follow perturbations of the tem- 
perature. Following DCHLC99 we write 

T— = Fiir — ^^cmit J (10) 

at 

where F is the thermal heat capacity. After substituting 
Eqs.(||,||j^), the definitions of Fi^,. and -Fcmitj and defin- 
ing the perturbation Tc{r,t) = T°{r){l + tl){r,t)r'^), where 
Tg (r) is the temperature belonging to the static flaring 
disk solution Eq.(H), one finds the following perturbation 
equation (see DCHLC99 for details): 



9-0 
'dt 



C 



dr 



(11) 



where t = {a^T^ /T)t. Eq.(pT|) is an advection equation. 
Apparently a perturbation moves inwards, and its ampli- 
tude damps, since is conserved in amplitude along the 
inward moving characteristic, and therefore the comoving 
amplitude of ijir"^ becomes smaller (DCHLC99). 

3.2. Hydrodynamic perturbations at thermal equilibrium 

When the cooling time of the disk is much shorter than 
the orbital time, one can assume that the disk is always 
in thermal equilibrium, and the disk height will react to 
changes in the irradiation flux in a hydrodynamic way, 
instead of being in vertical hydrostatic balance. In prin- 
ciple one should solve the time-dependent vertical hydro- 
dynamics of the disk (1-D hydrodynamics). But for the 
present analysis we simplify the picture, and assume that 



the density profile is always given by Eq.(^, but with 
a function of time. This assumption can be justified, since 
in the early (linear) stages of a possible instability one 
has dvz/dt oc z, which preserves the Gaussian shape of 
the density profile. By neglecting quadratic terms in the 
vertical velocity, one can then derive 



d^H„ 



kn 1 



dt^ MgWu 



(12) 



When we define the dimensionless time t = y/ GM^J R^t, 
a dimensionless form of the above equation is 



d^h{r,i) 
dP 



r^/^h(r,i) 



dh{r,i) 
dr 



1/4 



h{r, i) 



(13) 



with G as defined in Eq.(^ 

We now analyze the time-dependent behavior of linear 
perturbations of the power law solution (Eq.0) . We take 



Hr,i) 



7C 



1/7 



r^/^ [l + ^j{r,i)]. 



(14) 



Substituting this in Eq.(|l3[) yields (exactly): 



5V 



r3 1(1 + ip) 



1+-0 + 



7 &tp_ 



1/4 



(1 + V^) -(15) 



For tp <^ 1 this becomes 



(9V 



7 



di^ 8r3 i dr 



This equation has solutions of the form 
iP = Ar^e""* sin (uji - kr^ + (pa) , 
with the following dispersion relations: 

2 2 

uj — a , 
k = -(16/21) wcr. 



(16) 



(17) 



(18) 
(19) 



For inwards propagating waves (i.e. k/to < 0) one has 
CT > 0, and hence an exponentially growing modqj. The in- 
stability growth rate is inversely proportional to the wave- 
length of the perturbation, which indicates that the insta- 
bility will be dominated by the shortest wavelengths. It 
should be noted that the validity of the equations breaks 
down at the scale of the disk vertical height. For modes 
with larger k one should take radial radiative diffusion into 
account. This is, however, beyond the scope of the present 
set of equations. 



^ After submission the author became aware of work by 
E. Chiang (to be published) which reports on similar results. 
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4. Conditions for instability 



5. Discussion 



The instability only sets in under certain physical condi- 
tions. The most obvious conditions for this instability to 
occur are: 

1. Radiative cooling/heating time scale short compared 
to vertical sound crossing time, i.e. tdya 3> itherm- 

2. Thermal dust-gas coupling is strong enough to keep 
the gas cooling time shorter than the vertical sound 
crossing time. 

3. The equatorial temperature is dominated by irradia- 
tion, i.e. viscous dissipation by accretion is small in 
comparison. 

4. The disk is optically thick to starlight in radial (graz- 
ing) direction. 

5. The irradiation occurs through flaring; i.e. the reflected 
stellar flux from a scattering corona is weak in compar- 
ison to the direct interception of stellar light through 
flaring. 

We have not investigated whether disk models that obey 
all the above points, but are vertically optically thin in 
the near/mid infrared, are also unstable. 

Condition ^ is fulfilled for R ^ i?inst , where i?inst can 
be derived by equating the vertical sound crossing time 
with the thermal time: <dyn — ithorm- We have 



^■dyn 




(20) 



(21) 



2mu(TR.r| ' 

The surface density S(-R) (in units of g/cm^) is a free 
function of the model, as long as the disk remains optically 
thick to the incident stellar radiation. Equating tdyn — 
ithcrm yields the radius i?inst beyond which the instability 
may set in: 

13/3 



AU 



where L» is the luminosity of the star. 

The instability found here is an instability of the highly 
simplified analytic flaring disk models. Including more re- 
alistic physics could perhaps stabilize the disk. Some el- 
ements of realistic physics that should be considered in 
future work before one can claim that the instability is 
real are: 

1. 2-D axisymmetric hydrodynamics. 

2. Radial radiative transfer of stellar radiation instead of 
the inclination angle formula. 

3. 1-D vertical radiative transfer, because the top layers 
have a higher temperature (see CG97). 

4. Radial diffusion of radiation, or even fully 2-D radia- 
tive transfer, because radial exchange of energy might 
counter-act the shadowing effect which causes the in- 
stability. 



This Letter presents an analysis of the structure and sta- 
bility of irradiated non-accreting protostellar disks based 
on a highly simplified set of equations. In the context 
of these equations it is concluded that the flaring non- 
accreting reprocessing disk model is unstable when the 
vertical sound crossing time is larger than the heat- 
ing/cooling time scale. It should be noted, however, that 
detailed multidimensional numerical modeling is required 
to verify whether this instability is real or merely an arti- 
fact of the over-simplifled equations. 

The instability has a simple physical interpretation. 
Consider a flaring disk solution, and perturb it with an 
enlargement of the scale height at some point. The 'sunny 
side' of the hill will get overheated and expand verti- 
cally in a run-away fashion. The 'shadow side' will receive 
insufficient radiation and collapse. As the perturbations 
grow into the non-linear regime they start to cast a 'real' 
shadow {a < 0) over the disk at larger radii, thus com- 
pletely depriving this part of the disk of irradiation. At 
this point the validity of the simplified equations breaks 
down, and what happens in this non-linear regime is un- 
clear. Disk self irradiation (e.g. Bell 1999 ) might become 
important at this stage. 

The reason why we find an instability, while DCIILC99 
found only damping modes is because they considered only 
perturbations in the disk where the thermal time scale ex- 
ceeds the vertical sound crossing time scale. While for ac- 
tively accreting disks this condition is always satisfied, for 
passive disks the reverse may be true. This was recognized 
by DCHLC99, but an analysis of modes on a dynamic time 
scale in the outer regions of the protostellar disks was not 
carried out, hence their different conclusion. 
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